home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-07-11  |  4.8 KB  |  155 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  */
  17.  
  18. #ifndef lint
  19. static char sccsid[] = "@(#)log.c    5.2 (Berkeley) 4/29/88";
  20. #endif /* not lint */
  21.  
  22. /* LOG(X)
  23.  * RETURN THE LOGARITHM OF x 
  24.  * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
  25.  * CODED IN C BY K.C. NG, 1/19/85;
  26.  * REVISED BY K.C. NG on 2/7/85, 3/7/85, 3/24/85, 4/16/85.
  27.  *
  28.  * Required system supported functions:
  29.  *    scalb(x,n)
  30.  *    copysign(x,y)
  31.  *    logb(x)    
  32.  *    finite(x)
  33.  *
  34.  * Required kernel function:
  35.  *    log__L(z) 
  36.  *
  37.  * Method :
  38.  *    1. Argument Reduction: find k and f such that 
  39.  *            x = 2^k * (1+f), 
  40.  *       where  sqrt(2)/2 < 1+f < sqrt(2) .
  41.  *
  42.  *    2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  43.  *         = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  44.  *       log(1+f) is computed by
  45.  *
  46.  *                 log(1+f) = 2s + s*log__L(s*s)
  47.  *       where
  48.  *        log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
  49.  *
  50.  *       See log__L() for the values of the coefficients.
  51.  *
  52.  *    3. Finally,  log(x) = k*ln2 + log(1+f).  (Here n*ln2 will be stored
  53.  *       in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact
  54.  *       since the last 20 bits of ln2hi is 0.)
  55.  *
  56.  * Special cases:
  57.  *    log(x) is NaN with signal if x < 0 (including -INF) ; 
  58.  *    log(+INF) is +INF; log(0) is -INF with signal;
  59.  *    log(NaN) is that NaN with no signal.
  60.  *
  61.  * Accuracy:
  62.  *    log(x) returns the exact log(x) nearly rounded. In a test run with
  63.  *    1,536,000 random arguments on a VAX, the maximum observed error was
  64.  *    .826 ulps (units in the last place).
  65.  *
  66.  * Constants:
  67.  * The hexadecimal values are the intended ones for the following constants.
  68.  * The decimal values may be used, provided that the compiler will convert
  69.  * from decimal to binary accurately enough to produce the hexadecimal values
  70.  * shown.
  71.  */
  72.  
  73. #if defined(vax)||defined(tahoe)    /* VAX D format */
  74. #include <errno.h>
  75. #ifdef vax
  76. #define _0x(A,B)    0x/**/A/**/B
  77. #else    /* vax */
  78. #define _0x(A,B)    0x/**/B/**/A
  79. #endif    /* vax */
  80. /* static double */
  81. /* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
  82. /* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
  83. /* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
  84. static long     ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)};
  85. static long     ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)};
  86. static long     sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)};
  87. #define    ln2hi    (*(double*)ln2hix)
  88. #define    ln2lo    (*(double*)ln2lox)
  89. #define    sqrt2    (*(double*)sqrt2x)
  90. #else    /* defined(vax)||defined(tahoe)    */
  91. static double
  92. ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
  93. ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
  94. sqrt2  =  1.4142135623730951455E0     ; /*Hex  2^  0   *  1.6A09E667F3BCD */
  95. #endif    /* defined(vax)||defined(tahoe)    */
  96.  
  97. double log(x)
  98. double x;
  99. {
  100.     static double zero=0.0, negone= -1.0, half=1.0/2.0;
  101.     double logb(),scalb(),copysign(),log__L(),s,z,t;
  102.     int k,n,finite();
  103.  
  104. #if !defined(vax)&&!defined(tahoe)
  105.     if(x!=x) return(x);    /* x is NaN */
  106. #endif    /* !defined(vax)&&!defined(tahoe) */
  107.     if(finite(x)) {
  108.        if( x > zero ) {
  109.  
  110.        /* argument reduction */
  111.           k=logb(x);   x=scalb(x,-k);
  112.           if(k == -1022) /* subnormal no. */
  113.            {n=logb(x); x=scalb(x,-n); k+=n;} 
  114.           if(x >= sqrt2 ) {k += 1; x *= half;}
  115.           x += negone ;
  116.  
  117.        /* compute log(1+x)  */
  118.               s=x/(2+x); t=x*x*half;
  119.           z=k*ln2lo+s*(t+log__L(s*s));
  120.           x += (z - t) ;
  121.  
  122.           return(k*ln2hi+x);
  123.        }
  124.     /* end of if (x > zero) */
  125.  
  126.        else {
  127. #if defined(vax)||defined(tahoe)
  128.         extern double infnan();
  129.         if ( x == zero )
  130.             return (infnan(-ERANGE));    /* -INF */
  131.         else
  132.             return (infnan(EDOM));    /* NaN */
  133. #else    /* defined(vax)||defined(tahoe) */
  134.         /* zero argument, return -INF with signal */
  135.         if ( x == zero )
  136.             return( negone/zero );
  137.  
  138.         /* negative argument, return NaN with signal */
  139.         else 
  140.             return ( zero / zero );
  141. #endif    /* defined(vax)||defined(tahoe) */
  142.         }
  143.     }
  144.     /* end of if (finite(x)) */
  145.     /* NOTREACHED if defined(vax)||defined(tahoe) */
  146.  
  147.     /* log(-INF) is NaN with signal */
  148.     else if (x<0) 
  149.         return(zero/zero);      
  150.  
  151.     /* log(+INF) is +INF */
  152.     else return(x);      
  153.  
  154. }
  155.